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We present a realistic treatment of the hydrodynamic evolution of ultrarelativistic heavy ion 
collisions, based on the following features: initial conditions obtained from a flux tube approach, 
compatible with the string model and the color glass condensate picture; event-by-event procedure, 
taking into the account the highly irregular space structure of single events, being experimentally 
visible via so-called ridge structures in two-particle correlations; use of an efficient code for solv- 
ing the hydrodynamic equations in 3+1 dimensions, including the conservation of baryon number, 
strangeness, and electric charge; employment of a realistic equation-of-state, compatible with lattice 
gauge results; use of a complete hadron resonance table, making our calculations compatible with 
the results from statistical models; hadronic cascade procedure after an hadronization from the 
thermal matter at an early time. 



I. INTRODUCTION 

There seems to be little doubt that heavy ion collisions 
at RHIC energies produce matter which expands as an 
almost ideal fluid [lHj| ■ This observation is mainly based 
on the studies of azimuthal anisotropies, which can be 
explained on the basis of ideal hydrodynamics [^-Q. A 
big success of this approach was the correct description of 
the so-called mass splitting, which refers to quite different 
transverse momentum dependencies of the asymmetries 
for the different hadrons, depending on their masses. 

Another striking observation is the fact that particle 
production seems to be governed by statistical hadroniza- 
tion in the framework of an ideal resonance gas, with a 
hadronization temperatures Th close to 170 MeV (Tol " 
[l6| . which corresponds to the critical temperature of the 
(cross-over) transition between the resonance gas and the 
quark gluon plasma. Such a high temperature is in par- 
ticular necessary to accommodate the yields of heavy par- 
ticles like baryons and antibaryons. 

If we imposed statistical hadronization at Th ~ 170 
MeV in a hydrodynamical approach, we would get the 
correct particle ratios, but the baryon spectra would be 
too soft A later freeze-out at around 130 - 140 MeV, 
as in earlier calculations, gives better spectra, but too 
few baryons. A way out is to consider an early "chem- 
ical freeze-out" Teh ~ Th, and then force the particle 
yields to stay constant till the final "thermal freeze-out" 
Tth [13 • Although in this way one might be able to un- 
derstand particle yields and spectra, such an approach 
produces too much azimuthal asymmetry (expressed via 
the second Fourier coefficient 112) compared to the data, 
in particular at large rapidities. Here, it seems to help 



to replace the hydrodynamic treatment of the evolution 
between Tch and Tth by a hadronic cascade [T8| - |2]| . So 
this second phase seems to be significantly non-thermal. 

The calculations of (Tsl . IT9| manage to reproduce both 
particle yields and transverse momentum spectra of pi- 
ous, kaons, and protons within 30%, for pt values below 
1.5 GeV/c. The net baryon yield cannot be reproduced, 
since the calculations are done for zero baryon chemical 
potential, another systematic problem is due to a rela- 
tively small hadron set. A bigger hadron set will produce 
essentially more pious and will thus reduce for example 
the pion / kaon ratio. 

Most calculations are still done using an unrealistic 
equation-of-state with a first order transition, based on 
ideal gases of quarks & gluons and hadrons. As shown 
later, it actually makes a big difference using a realistic 
equation-of-state, which is for /is = compatible with 
lattice results. 

Also important is an explicit treatment of individual 
events rather than taking smooth initial conditions rep- 
resenting many events. This has been pioneered by Sphe- 
rio calculations (2342^ , based on Nexus initial conditions 
[25l . [26| . An event-by-event treatment will affect all ob- 
servables like spectra and elliptical flow, and it is abso- 
lutely essential for rapidity-angle correlations (ridge ef- 
fect). 

Although Nexus reproduces qualitatively the essential 
features of a realistic event- by- event initial condition, it 
should be noted that the model has been developed ten 
years ago, before the RHIC era. So we will base our 
discussions in this paper on the Nexus successor EPOS, 
which contains many upgrades, related to the question 
of the interplay between soft and hard physics, high 
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partoii density effects and saturation, the role of pro- 
jectile and target remnants, and so on. The param- 
eters have been optimized by comparing to all possi- 
ble accelerator data concerning proton-proton (or more 
generally hadron-proton) and proton-nucleus (deuteron- 
nucleus) collisions. EPOS seems to be the only model 
compatible with yields, spectra, and double differential 
spectra of identified particles from NA49 [13 • EPOS 
seems as well to be the only interaction model compatible 
with cosmic ray data for air shower simulations |28| . All 
this just to say that we consider the elementary EPOS 
model for pp scattering as a very solid basis for general- 
izations towards heavy ion applications. 

In this paper, we present a realistic treatment of the 
hydrodynamic evolution of ultrarelativistic heavy ion col- 
lisions, based on the following features: 

• initial conditions obtained from a flux tube ap- 
proach (EPOS), compatible with the string model 
used since many years for elementary collisions 
(electron-positron, proton proton), and the color 
glass condensate picture; 

• consideration of the possibility to have a (moderate) 
initial collective transverse flow; 

• event-by-event procedure, taking into the account 
the highly irregular space structure of single events, 
being experimentally visible via so-called ridge 
structures in two-particle correlations; 

• core-corona separation, considering the fact that 
only a part of the matter thermalizes; 

• use of an efficient code for solving the hydrodynamic 
equations in 3+1 dimensions, including the conser- 
vation of baryon number, strangeness, and electric 
charge; 

• employment of a realistic equation-of-state, com- 
patible with lattice gauge results - with a cross-over 
transition from the hadronic to the plasma phase; 

• use of a complete hadron resonance table, making 
our calculations compatible with the results from 
statistical models; 

• hadronic cascade procedure after hadronization 
from the thermal system at an early stage. 

All the above mentioned features are not new, what 
is new is the attempt to put all these elements into a 
single approach, bringing together topics like statistical 
hadronization, flow features, saturation, the string model, 
and so on, which are often discussed independently. For 
any quantitative analysis of heavy ion results we have to 
admit that there is just one common mechanism, which 
accounts for the whole soft physics. We therefore test 



our approach by comparing to all essential observables in 
Au-Au scatterings at RHIC. 

There is quite some activity concerning viscous effects 
[29143^ . but this aspect will not be addressed in the 
present paper. Here, we want to develop a realistic de- 
scription based on ideal hydrodynamics, and see how far 
one can get. As we will see later, some of the features 
attributed to viscosity may be explained within ideal hy- 
drodynamics, in a realistic formulation. 

Although the model is very complex, the physical pic- 
ture which emerges is very clear , since the different "fea- 
tures" of our approach affect different observables in a 
very transparent way. A gold-gold collision at 200 GeV 
will typically create after less than one fm/c thermalized 
quark/gluon matter, concentrated in several longitudinal 
sub-flux-tube with energy density maxima of well beyond 
50 GeV/fm"^. Flux-tube structure essentially means that 
the complicated bumpy transverse structure of a given 
event is (up to a factor) translational invariant. During 
the evolution, translational invariant flows develop, which 
finally show up as rapidity-angle correlations. This is un- 
avoidable in such an approach with irregular flux tubes. 

In fig. [TJ we sketch the flux-tube picture.. The Ion- 




Figure 1: Macroscopic flux tubes (three in this example), 
made out of many individual ones, of variable length. 

gitudinal direction is along the z-axis, the coordinates x 
and y represent the transverse plane. A "macroscopic" 
flux tube is a longitudinal structure of high energy den- 
sity, almost translational invariant despite an irregular 
form in transverse direction. Such a flux tube is made 
of many individual elementary flux-tubes or strings, each 
on having a small diameter (of 0.2 to 0.3 fm). The el- 
ementary flux tubes are actually short, the momentum 
fraction of the string ends are distributed roughly as 1/x. 
The macroscopic flux tubes represent nevertheless long 
structures, simply due to the fact that many short el- 
ementary flux tubes are located at transverse positions 
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corresponding to the positions of nucleon-nucleon scat- 
terings. And these simply happen to be more or less 
frequent in certain transverse areas, as indicated in the 
figure by the three clusters of interaction positions (dots 
in the x ^ y plane). 

This flux tube approach is just a continuation of 30 
years of very successful applications of the string ap- 
proach to particle production in collisions of high en- 
ergy particles |35l437| . in particular in connection with 
the parton model. Here, the relativistic string is a phe- 
nomenological tool to deal with the longitudinal char- 
acter of the final state partonic system. An important 
issue at high energies is the appearance of so-called non- 
linear effects, which means that the simple linear par- 
ton evolution is no longer valid, gluon ladders may fuse 
or split. More recently, a classical treatment has been 
proposed, called Color Glass Condensate (CGC), having 
the advantage that the framework can be derived from 
first principles (38l - [43j . Comparing a conventional string 
model like EPOS and the CGC picture: they describe 
the same physics, although the technical implementation 
is of course different. All realistic string model implemen- 
tations have nowadays to deal with screening and satu- 
ration, and EPOS is not an exception. Without screen- 
ing, proton-proton cross sections and multiplicities will 
explode at high energies. We will discuss later in more 
detail about the question of CGC initial conditions for hy- 
drodynamical evolutions compared to conventional ones. 
To give a short answer: this question is irrelevant when 
it comes to event by event treatment. 

Starting from the flux-tube initial condition, the sys- 
tem expands very rapidly, thanks to the realistic cross- 
over equation-of-state, flow (also elliptical one) develops 
earlier compared to the case a strong first order equation- 
of-state as in (Tsl . [l9| , temperatures corresponding to the 
cross-over (around 170 MeV) are reached in less than 
10 fm/c. The system hadronizes in the cross-over re- 



gion, where here "hadronization" is meant to be the end 
of the completely thermal phase: matter is transformed 
into hadrons. We stop the hydrodynamical evolution at 
this point, but particles are not yet free. Our favorite 
hadronization temperature is 166 MeV, shown as the dot- 
ted line in fig. [51 which is indeed right in the transition 
region, where the energy density varies strongly with tem- 
perature. At this point we employ statistical hadroniza- 
tion, which should be understood as hadronization of the 
quark-gluon plasma state into a hadronic system, at an 
early stage, not the decay of a resonance gas in equilib- 
rium. 

After this hadronization -although no longer thermal- 
the system still interacts via hadronic scatterings, still 
building up (elliptical) flow, but much less compared to 
an idealized thermal resonance gas evolution, which does 
not exist in reality. 

Despite the non-equilibrium behavior in the finale stage 
of the collision, our sophisticated procedure gives parti- 
cle yields close to what has been predicted in statistical 
models, see fig. [3] This is because the final hadronic cas- 
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Figure 2: The energy density over as a function of the 
temperature T. The dotted line indicates the "hadronization 
temperature", i.e. end of the thermal phase, when "matter" is 
transformed into hadrons. 



Figure 3: Particle ratios (hadron yields to tv^ yields) from our 
model calculations (thick horizontal line) as compared to the 
statistical model [lo| (thin horizontal line), and to data (43 - 14^ 1 
(points). 

cade does not change particle yields too much (with some 
exceptions to be discussed later), but it affects slopes and 
-as mentioned- azimuthal asymmetry observables. 

In the following, we will present the details of our re- 
alistic approach to the hydrodynamic evolution in heavy 
ion collisions, with a subsequent attempt to understand 
and interpret all soft heavy ion data from Au-Au at 200 
GeV. The predictive power of the presented approach is 
enormous. The basic EPOS approach, which fixes the 
fiux tube initial conditions, has quite a number of pa- 
rameters determining soft Pomeron properties, the per- 
turbative QCD treatment (cutoffs), the string dynamics, 
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screening and saturation effects, the projectile and target 
remnant properties. All these unknowns are fixed by in- 
vestigating electron-positron, proton-proton, and proton- 
nucleus scattering from SPS via RHIC to Tevatron ener- 
gies, for all observables where data are available. This 
huge amount of elementary data lets very little freedom 
concerning heavy ion collisions. 



II. ELEMENTARY FLUX TUBES AND 
NON-LINEAR EVOLUTION 

Nucleus-nucleus scattering - even proton-proton - 
amounts to many elementary collisions happening in par- 
allel. Such an elementary scattering is the so-called "par- 
ton ladder" , see fig. HI also referred to as cut Pomeron, 
see appendix 1X1 and |47| . A parton ladder represents par- 
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Figure 4: Elementary interaction in the EPOS model. 

ton evolutions from the projectile and the target side to- 
wards the center (small x) . The evolution is governed by 
an evolution equation, in the simplest case according to 
DGLAP. In the following we will refer to these partons as 
"ladder partons", to be distinguished from "spectator par- 
tons" to be discussed later. It has been realized a long 
time ago that such a parton ladder may be considered 
as a quasi-longitudinal color field, a so-called "flux tube", 
conveniently treated as a relativistic string. The inter- 
mediate gluons are treated as kink singularities in the 
language of relativistic strings, providing a transversely 
moving portion of the object. This flux tube decays via 
the production of quark-antiquark pairs, creating in this 
way fragments - which are identifled with hadrons. Such 
a picture is also in qualitative agreement with recent de- 
velopments concerning the Color Glass Condensate, as 
discussed earlier. 



A consistent quantum mechanical treatment of the 
multiple scattering is quite involved, in particular when 
the energy sharing between the parallel scatterings is 
taken into account. For a detailed discussion we refer 
to [1^. Based on cutting rule techniques, one obtains 
partial cross sections for exclusive event classes, which 
are then simulated with the help of Markov chain tech- 
niques. 

Important in particular at moderate energies (RHIC): 
our "parton ladder" is meant to contain two parts |25| : 
the hard one, as discussed above (following an evolution 
equation), and a soft one, which is a purely phenomeno- 
logical object, parametrized in Regge pole fashion, see 
appendix. The soft part essentially compensates for the 
infrared cutoffs, which have to be employed in the per- 
turbative calculations. 

At high energies, one needs to worry about non-linear 
effects, because the gluon densities get so high that gluon 
fusion becomes important. Nonlinear effects could be 
taken into account in the framework of the CGC (ssl - 
|43| . Here , we adopt a phenomenological approach, which 
grasps the main features of these non-linear phenomena 
and still remains technically doable (we should nor forget 
that we finally have to deal with complications due to 
multiple scatterings, as discussed earlier). 

Our phenomenological treatment is based on the fact 
that there are two types of nonlinear effects jl^]: a sim- 
ple elastic rescattering of a ladder parton on a projectile 
or target nucleon (elastic ladder splitting) . or an inelastic 
rescattering (inelastic ladder splitting), see figs. [5l|6l The 
elastic process provides screening, therefore a reduction of 
total and inelastic cross sections. The importance of this 
effect should first increase with mass number (in case of 
nuclei being involved), but finally saturate. The inelas- 
tic process will affect particle production, in particular 
transverse momentum spectra, strange over non-strange 
particle ratios, etc. Both, elastic and inelastic rescat- 
tering must be taken into account in order to obtain a 
realistic picture. 

nucleons 




ladder partons 

Figure 5: Elastic "rescattering" of a ladder parton. 

To include the effects of elastic rescattering, we first 
parametrize a parton ladder (to be more precise: the 
imaginary part of the corresponding amplitude in impact 
parameter space) computed on the basis of DGLAP. We 
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Figure 6: Inelastic "rescattering" of a ladder parton. 

obtain an excellent fit of the form 



(1) 



where and x~ are the momentum fractions of the 
"first" ladder partons on respectively projectile and target 
side (which initiate the parton evolutions). The parame- 
ters a and /3 depend on the cms energy of the hadron- 
hadron collision. To mimic the reduction of the increase 
of the expressions a{x'^x~)^ with energy, we simply re- 
place them by 



a{x+)'^+'^x-) 



(2) 



where the values of the positive numbers Sp/x will in- 
crease with the nuclear mass number and log s. This addi- 
tional exponent has very important consequences: it will 
reduce substantially the increase of both cross sections 
and multiplicity with the energy, having thus a similar 
effect as introducing a saturation scale. 

The inelastic rescatterings (ladder splittings, looking 
from insider to outside) amount to providing several lad- 
ders close to the projectile (or target) side, which are close 
to each other in space. They cannot be considered as in- 
dependent color fields (strings), we should rather think of 
a common color field built from several partons ladders. 
We treat this object via an enhancement of remnant ex- 
citations, the latter ones to be discussed in the following. 

So far we just considered two interacting partons, one 
from the projectile and one from the target. These 
partons leave behind a projectile and target remnant, 
colored, so it is more complicated than simply projec- 
tile/target deceleration. One may simply consider the 
remnants to be diquarks, providing a string end, but this 
simple picture seems to be excluded from strange an- 
tibaryon results at the SPS We therefore adopt the 
following picture: not only a quark, but a two-fold object 
takes directly part in the interaction, namely a quark- 
antiquark or a quark-diquark pair, leaving behind a col- 
orless remnant, which is, however, in general excited (off- 
shell). If the first ladder parton is a gluon or a seaquark, 
we assume that there is an intermediate object between 
this gluon and the projectile (target), referred to as soft 
Pomeron. And the "initiator" of the latter one is again 
the above-mentioned two-fold object. 



So we have finally three "objects", all of them being 
white: the two off-shell remnants, and the parton ladder 
in between. Whereas the remnants contribute mainly 
to particle production in the fragmentation regions, the 
ladders contribute preferentially at central rapidities. 

We showed in ref. ji^ that this "three object picture" 
can solve the "multi-strange baryon problem" of ref. [4^ . 
In addition, we assembled all available data on particle 
production in pp and pA collisions between 100 GeV (lab) 
up to Tevatron, in order to test our approach. Large ra- 
pidity (fragmentation region) data are mainly accessible 
at lower energies, but we believe that the remnant prop- 
erties do not change much with energy, apart of the fact 
that projectile and target fragmentation regions are more 
or less separated in rapidity. But even at RHIC, there are 
remnant contribution at rapidity zero, for example the 
baryon/antibaryon ratios are significantly different from 
unity, in agreement with our remnant implementation. 
So even central rapidity RHIC data allow to confirm our 
remnant picture. 



III. FLUX TUBES, JETS, AND CORE-CORONA 
SEPARATION 



We will identify parton ladders with elementary flux 
tubes, the latter ones treated as classical strings. The rel- 
ativistic classical string picture is very attractive, because 
its dynamics (Lagrangian) is essentially derived from gen- 
eral principles as covariance and gauge invariance (the dy- 
namics should not depend on a particular string surface 
parametrization). We use the simplest possible string: 
a two-dimensional surfaces X{a, (3) in 3+1 dimensional 
space-time, with piecewise constant initial conditions. 



r)X 

V{a)^—{a,P 



0) = Vfe, in [Q;fc,afc+i], (3) 



referred to as kinky strings. The dynamics is governed 
by the Nambu-Goto string action (50l - [5^ (see also (36|). 
Our string is characterized by a sequence of intervals 
[afc,afc+i], and the corresponding velocities Vfe. Such 
an interval with its constant value of V is referred to 
as "kink". Now we are in a position to map partons onto 
strings: we identify the ladder partons with the kinks of 
a kinky string, such that the length of the a-interval is 
given by the parton energies Ej^, and the kink velocities 
are just the parton velocities, p'^/E^. The string evolu- 
tion is then completely given by these initial conditions, 
expressed in terms of parton momenta. The string sur- 
face is given as 



X{a,P) = Xo + 



1 



a+i3 



a-[j 



(4) 



Let us considers a string at a given proper time tq. In 
fig. [71 the thick line of the form of a hyperbola repre- 
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Figure 7: Flux tube (string) at given proper time. Tiie picture 
is schematic in the sense that the string extends well into 
the transverse dimension, correctly taken into account in the 
calculations. The quantity X is a four-vector! 

sents schematically the intersection of the string surface 
X{a,(3) with the hypersurface corresponding to constant 
proper time: r = tq. We show only a simplified picture in 
z — t space, whereas in reality (and in our calculations) all 
three space dimensions are important, due to the trans- 
verse motion of the kinks: the string at constant proper 
time is a one-dimensional manifold in the full 3+1 dimen- 
sional space-time. In fig. \8\ we sketch the space compo- 
nents of this object: the string in M'^ space is a mainly 
longitudinal object (here parallel to the z-axis) but due to 
the kinks there are string pieces moving transversely (in 
y-direction in the picture). But despite these kinks, most 
of the string carries only little transverse momentum! 




Figure 8: Flux tube with transverse kink in IR"^ space. The 
kink leads to transversely moving string regions (transverse 
arrow) . 

In case of elementary reactions like electron-positron 
annihilation or proton proton scattering (at moderately 
relativistic energies), hadron production is realized via 
string breaking, such that string fragments are identified 
with hadrons. Here, we employ the so-called area law hy- 
pothesis [l^llJl (see also Hg): the string breaks via q — q 
ov qq — qq production within an infinitesimal area dA on 
its surface with a probability which is proportional to this 
area, dP = ps dA, where pb is the fundamental parame- 
ter of the procedure. It should be noted that despite the 
very complicated structure of the string surface X(a,/3) 
in 3+1 space-time, the breaking procedure following the 
area law can be done rigorously, using the so-called band- 
method [1^ [5^. The flavor dependence of the q — q 



ov qq — qq string breaking is given by the probabilities 
exp(— Trm^/ft), with ruq being the quark masses and k 
the string tension. After breaking, the string pieces close 
to a kink constitute the jets of hadrons (arrows in fig. [5]), 
whose direction is mainly determined by the kink-gluon. 




Figure 9: Broken flux tube with transverse kink in IR"^ space. 
The string segments close to the kink giving rise to trans- 
versely moving hadrons, constituting a jet (arrows). 

When it comes to heavy ion collisions or very high en- 
ergy proton-proton scattering, the procedure has to be 
modified, since the density of strings will be so high that 
they cannot possibly decay independently (56| . For tech- 
nical reasons, we split each string into a sequence of string 
segments, corresponding to widths 5a and Sp in the string 
parameter space (see fig. [TUl One distinguishes between 




" ^ X(a+5a,|3-H6P) 
X(a,p) 

Figure 10: String segment at given proper time. The picture 
is schematic in the sense that the string extends well into 
the transverse dimension, correctly taken into account in the 
calculations. 

string segments in dense areas (more than some critical 
density po of segments per unit volume), from those in 
low density areas. The high density areas are referred to 
as core, the low density areas as corona jlBl- String seg- 
ments with large transverse momentum (close to a kink) 
are excluded from the core. At this stage, we do not 
consider energy loss of these kink partons, we will inves- 
tigate this in a later publication. Also excluded from the 
core are remnant baryons. Simple imp lementations of the 
core-corona idea can be found in |57l [58|. 

Let us consider the core part. Based on the four- 
momenta of infinitesimal string segments, we compute 
the energy momentum tensor and the flavor flow vector 
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at some position x (at t = tq) as [2^ 

where q g u,d, s represents the net flavor content of the 
string segments, and 

are the four-momenta of the segments. The function gr is a 
Gaussian smoothing kernel with a transverse width cr± = 
0.25 fm. The Lorentz transformation into the comoving 
frame gives 

A" ^A^ .T^"^ = TZn, (8) 

where we define the comoving frame such that the first 
column of Tcom is of the form (e, 0, 0, 0)^. This provides 
an equation for the energy density e in the comoving 
frame, and the flow velocity components : 

3 

k=l 

which may be solved iteratively (s^ . 

3 

k=l 

^(n)i ^ J_(-j,*0_yifc^(«-l)fcN 

£(") ^ / \ J 

The flavor density is then calculated as 

fg = NgU, (13) 

with u being the flow four- velocity. 

From the above procedure, we get event- by- event fluc- 
tuations of the collective transverse velocities, but these 
flows are very small. However, several authors [60 65] 
discussed recently the possibility of having already an 
initial collective velocity. We consider such a possibil- 
ity by adding to our transverse velocities v^/y (r, 4>) the 
following terms: 

Avx{r,(f)) = min(0.4, VQr/ro) \/l + e coscf), (14) 
Avy{r,(f>) = min(0.4, Vor/ro) — e sin^, (15) 

with 



and 

p-4v/(x2 + y2)/2, e^{y^^x^)/{y'+x'). (17) 

Such an initial collective transverse flow seems to be not 
really essential for reproducing the data, however, a value 
vq = 0.25 gives a slight improvement of the transverse 
momentum spectra, compared to wq = 0. So we use the 
former value as default. 



IV. HYDRODYNAMIC EVOLUTION, 
REALISTIC EQUATION-OF-STATE 

Having fixed the initial conditions, the core evolves ac- 
cording to the equations of ideal hydrodynamics, namely 
the local energy-momentum conservation 

a^r^'^ = o, zy = o,...,3, (18) 

and the conservation of net charges, 

dNj^ = 0, k = B,S,Q, (19) 

with B, S, and Q referring to respectively baryon num- 
ber, strangeness, and electric charge. In this paper we 
treat ideal hydrodynamic, so we use the decomposition 

T"'' = {e+p)u''u'' -pg"" , (20) 

K = UkU^, (21) 

where u is the four- velocity of the local rest frame. Solv- 
ing the equations, as discussed in the appendix, provides 
the evolution of the space-time dependence of the macro- 
scopic quantities energy density e{x), collective flow ve- 
locity v{x), and the net flavor densities nk{x). Here, the 
crucial ingredient is the equation of state, which closes 
the set of equations by providing the e-dependence of the 
pressure p. The equation-of-state should fulflU the fol- 
lowing requirements: 

• flavor conservation, using chemical potentials fiB, 
Ms, mq; 

• compatibility with lattice gauge results in case of 
Ms =M5 =MQ 0. 

The starting point for constructing this "realistic" 
equation-of-state is the pressure pn of a resonance gas, 
and the pressure pQ of an ideal quark gluon plasma, in- 
cluding bag pressure. Be Tc the temperature where pn 
and pq cross. The correct pressure is assumed to be of 
the form 



ro = p\/ 1 — e cos 2(j), 



(16) 



p = PQ + X{PH - Pq), (22) 
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where the temperature dependence of A is given as 
A = exp 



^ '^"]e{T-T,) + eiT,-T), (23) 



S J 



with 



S = So exp {-{flB/^^c)^) (^1 + 



T-Tc 



(24) 



From the pressure one obtains the entropy density S as 
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S =^ = Sq^\{Sh - Sq) + ^{vh-vq), (25) 



and the flavor densities as 

d-p 



dX 



= 'dl^^^'^Q^^^''H-nQ) + ^^{pH-PQ). (26) 



The energy density is finally given as 



(27) 



or 



s = eQ + X{eH~eQ)+( r|A ) (p^-pQ). (28) 



Our favorite equation-of-state, referred to as "X3F", is 
obtained for Sq = 0.15, which reproduces lattice gauge 
results for fis —fJ-s —fJ-Q = Oj shown in figs. 1111 and 



1- 22.5 I 




QIF 



^TT* • X3F 



0.5 1 1.5 



2 2.5 
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Figure 11: Energy density versus temperature, for our 
equation-of-state X3F (full line), compared to lattice data [6^ 
(points), and some other EoS choices, see text. 



Figure 12: Pressure versus temperature, for our equation-of- 
state X3F (full line), compared to lattice data [6y| (points), 
and some other EoS choices, see text. 



The symbol X3F stands for "cross-over" and "3 flavor 
conservation". Also shown in the figures is the EoS QIF, 
referring to a simple first order equation-of-state, with 
baryon number conservation, which we will use as a ref- 
erence to compare with. Many current calculations are 
still based on this simple choice, as for example the one 
in [Tsl. [l9|. shown as dotted lines in figs. [TTjandlT^ 

When the evolution reaches the hadronization hyper- 
surface, defined by a given temperature Th, we switch 
from "matter" description to particles, using the Cooper- 
Frye description. Particles may still interact, as dis- 
cussed below, so hadronization here means an interme- 
diate stage, particles are not yet free streaming, but they 
are not thermalized any more. The hadronization pro- 
cedure is described in detail in the appendix. After 
the "intermediate" hadronization, the particles at their 
hadronization positions (on the corresponding hypersur- 
face) are fed into the hadronic cascade model UrQMD 
[stI . [68| , performing hadronic interaction until the system 
is so dilute that no interaction occur any more. The "fi- 
nal" freeze out position of the particles is the last interac- 
tion point of the cascade process, or the hydro hadroniza- 
tion position, if no hadronic interactions occurs. 



V. ON THE IMPORTANCE OF AN 
EVENT-BY-EVENT TREATMENT 

A remarkable feature of an event-by-event treatment of 
the hydrodynamical evolution based on random fiux tube 
initial conditions is the appearance of a so-called ridge- 
structure, found in Spherio calulations based on Nexus 
initial conditions (6^, [73 . We expect to observe a similar 
structure doing an event- by- event hydrodynamical evo- 
lution based on fiux-tube initial conditions from EPOS. 
The result is shown in fig. 1131 where we plot the di- 
hadron correlation dN/dArjdAcI), with Arj and A(j> 
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Figure 13: Dihadron Arj — A(j> correlation in a central Au-Au collision at 200 GeV, as obtained from an event-by-event 
treatment of the hydrodynamical evolution based on random flux tube initial conditions. Trigger particles have transverse 
momenta between 3 and 4 GeV/c, and associated particles have transverse momenta between 2 GeV/c and the pt of the trigger. 



lergy density [GeV/fm3] (eta_s= 0.0 , tau= 0.6fm) C 1 





Figure 14: Initial energy density in a central Au-Au collision 
at 200 GeV, at a space-time rapidity 77s — 0. 



Figure 15: Initial energy density in a central Au-Au collision 
at 200 GeV, at a space-time rapidity ris = 1-5. 



being respectively the difference in pseudorapidity and 
azimuthal angle of a pair of particles. Here, we consider 
trigger particles with transverse momenta between 3 and 
4 GeV/c, and associated particles with transverse mo- 
menta between 2 GeV/c and the pt of the trigger, in cen- 
tral Au-Au collisions at 200 GeV. Our ridge is very similar 
to the structure observed by the STAR collaboration (7l| . 
In the following we will discuss a particular event. 



which can, however, be considered as a typical example, 
with similar observations being true for randomly cho- 
sen events. Important for understanding the strong A77 - 
A(f) correlation is the observation, that the initial energy 
density has a very bumpy structure as a function of the 
transverse coordinates x and y. However, this irregular 
structure is the same at different longitudinal positions. 
This can be clearly seen in figs. [T3] and [151 where we 
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energy density [GeV/fm3] (eta_5= 0.0 , tau= 2.6fm) C 1 




x[fni] 



Figure 16: Energy density at a proper time r — 2.6fm/c, at 
a space-time rapidity 77s — 0. 




Figure 17: Energy density at a proper time r — 2.6fm/c, at 
a space-time rapidity 77^ — 1.5. 



show for a given event the energy density distributions 
in the transverse planes at different space-time rapidi- 
ties, namely rjg ~ and rjs = 1-5: we observe almost the 
same structure. For different events, the details of the 
bumpy structures change, but we always find an approxi- 
mate "translation invariance": the distributions of energy 
density in the transverse planes vary only little with the 
longitudinal variable 77^. It should be noted that the col- 
ored areas represent only the interior of the hadronization 
surface, the outside regions are white. Hadronization is 
meant to be an intermediate step, before the hadronic 
cascade. An approximate translational invariance is also 
observed when we go to larger values of rjs , so for exam- 
ple when we compare the energy density at 77^ = 1.5 with 



rad velocity [% of c] {eta_s= 0.0 , tau= 2.6fm) C 1 




Figure 18: Radial flow velocity at a proper time r = 2.6fm/c, 
at a space-time rapidity i]s = 0. 



rad veiocity [% of c] {eta_s= 1 .5 , tau= 2.6fm) C 1 




Figure 19: Radial flow velocity at a proper time r = 2.6fm/c, 
at a space-time rapidity i]s = 1.5. 



the one at 77^ = 3.0: the form of the energy distributions 
is similar, however, the magnitude at large 77s is smaller. 

Considering later times, we see in figs. (TH] to [121 that 
the approximate translational invariance is conserved, for 
both energy densities and radial fiow velocities. It is re- 
markable (and again true in general, for arbitrary events) 
that the energy distribution in the transverse plane is 
much smoother than initially, the distribution looks more 
homogeneous. Very important for the following discus- 
sion is the fiow pattern, seen in figs. [TS]and[Tni for 77,, = 
and ?7s = 1.5 : the radial fiow is as expected largest in 
the outer regions. Closer inspection of the outside ring of 
large radial flows reveals an irregular atoll-like structure: 
there are well pronounced peaks of large fiow over the 



Event-by-Event Simulation of the Three-Dimensional Hydrodynamic Evolution 



11 



energy density [GeV/fm3] (eta_5= 0.0 , tau= 4.6fm) C 1 




Figure 20: Energy density at a proper time r — 4.6fm/c, at 
a space-time rapidity 77s — 0. 



energy density [GeV/fm3] (eta_s= 1 .5 , tau= 4.6f m) C 1 




x[fm] 



Figure 21: Energy density at a proper time r — 4.6fm/c, at 
a space-time rapidity 77s — 1.5. 



background ring. At even later times, as seen in figs. [50] 
todSl the outer surfaces get irregular, due to the irregular 
flows discussed above, again with well identified peaks of 
large radial flows. 

The well isolated peaks of the radial flow velocities have 
two important properties: they sit close to the hadroniza- 
tion surface, and they sit at the same azimuthal angle, 
when comparing different longitudinal positions rjg. As 
a consequence, particles emitted from different longitudi- 
nal positions get the same transverse boost , when their 
emission points correspond to the azimuthal angle of a 
common flow peak position. And since longitudinal co- 
ordinate and (pseudo)rapidity are correlated, one obtains 
finally a strong A77 - A(j> correlation. 



rad velocity [% of c] {eta_s= 0.0 , tau= 4.6fm) C 1 




Figure 22: Radial flow velocity at a proper time r = 4.6fm/c, 
at a space-time rapidity i]s = 0. 



rad velocity [% of c] {eta_s= 1 .5 , tau= 4.6fm) C 1 




Figure 23: Radial flow velocity at a proper time r = 4.6fm/c, 
at a space-time rapidity i]s = 1.5. 



In fig. [Ml we summarize the above discussion: the 
flux tube initial conditions provide a bumpy structure of 
the energy density in the transverse plane, which shows, 
however, an approximate translational invariance (simi- 
lar behavior at different longitudinal coordinates). Solv- 
ing the hydrodynamic equations preserves this invariance, 
leading in the further evolution to an invariance of the 
transverse flow velocities. These identical flow patterns 
at different longitudinal positions lead to the fact that 
particles produced at different values of rjs profit from 
the same collective push, when they are emitted at an 
azimuthal angle corresponding to a flow maximum (indi- 
cated by the arrows in the figure). 

Finally we have to address the question, why we have a 
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Figure 24: Schematic view of the translational invariance of 
the initial energy density (a), leading to a corresponding in- 
variance of the transverse flow. We use the term "invariance" 
in the sense of a similarity transform: same shape, but differ- 
ent magnitude. The magnitude of the energy density at large 
rjs is of course smaller than the one at rjs — 0. 

irregular transverse structure with an approximate trans- 
lational invariance. The basic structure of EPOS is such 
that each individual nucleon-nucleon collision results in 
a projectile and target remnant, and two or more ele- 
mentary flux tubes (strings). The higher the energy the 
bigger the number of strings. Most of the energy of the 
reaction is carried by the remnants, the flux tubes cover 
only a limited range in rapidity, but their "lengths" (in 
rapidity) vary enormously. Nevertheless we obtain a very 
smooth variation of the energy density with the longi- 
tudinal coordinate rjs. This is due to the fact that the 



transverse positions of a string is given by the position of 
the nucleon pair, who's interaction gave rise the the for- 
mation of the flux tube. These "pair positions" fluctuate 
considerably, event- by- event, and one obtains typically a 
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Figure 25: Projection of the positions of nucleon-nucleon scat- 
tering to the transverse {x,y) plane, from a simulation of a 
semi-peripheral (6 = 8fm/c) Au-Au event at 200 GeV. 




Figure 26: Schematic view of the projection of the positions 
of nucleon-nucleon scattering to the transverse {x,y) plane, 
which defines "possible transverse positions" of the flux tubes, 
indicated by the thin lines. The actual flux tubes fluctuate 
concerning their longitudinal positions; a possible realization 
is shown by the thick lines. 

situation as shown in fig. where we plot the projection 
to the transverse plane of the positions of the interaction 
nucleon-nucleon pairs. The two circles representing two 
hard sphere nuclei is only added to guide the eye, for 
the calculations we use of course a realistic nuclear den- 
sity. Clearly visible in the figure is the inhomogeneous 
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structure: there are areas with a high density of interac- 
tion points, and areas which are less populated. These 
transverse positions of interacting pairs define also the 
corresponding positions of the flux tubes associated to 
the pairs. In fig. 1261 we present a schematic view of this 
situation: on the left we plot the pair positions projected 
to the transverse plane (dots) . From each dot we draw a 
line parallel to the z-axis, representing a possible location 
of a fiux tube. The flux tubes have variable longitudinal 
lengths, they do not cover the full possible length between 
projectile and target, but only a portion, as indicated by 
the thick horizontal lines in the figure. But even then, the 
transverse structure (minima and maxima of the energy 
density) is to a large extend determined by the density of 
nucleon-nucleon pairs. 



VI. ELLIPTICAL FLOW 

Important information about the space-time evolution 
of the system is provided by the study of the azimuthal 
distribution of particle production. One usually expands 



dn 

— (X I + 2v2 cos 20 



(29) 



where a non-zero coefficient V2 is referred to as elliptical 
flow [t^- It is usually claimed that the elliptical flow is 
proportional to the initial space eccentricity 



(30) 



(y2 -I- ■ 

We therefore plot in fig. [^7] the ratio of V2 over ec- 




50 100 150 200 250 300 350 400 
Npart 

Figure 27: Centrality dependence of the ratio of V2 over ec- 
centricity. Points are data (7§ |. the different curves refer to 
the full calculation - hydro & cascade (full line), only elastic 
hadronic scatterings (dotted), and no hadronic cascade at all 
(dashed). The thin solid line -above all others- refers to the 
hydrodynamic calculation till final freeze-out at 130 MeV. 

centricity. The points are data; the full line is the full 
calculation: hydrodynamical evolution with subsequent 
hadronic cascade, from flux tube initial conditions, in 



event- by- event treatment. The dotted line refers to a sim- 
plified hadronic cascade, allowing only elastic scatterings, 
the dashed line is the calculation without hadronic cas- 
cade. In all cases, hadronization from the thermal phase 
occurs at Th = 166 MeV. We also show as thin solid line 
the hydrodynamic calculation till final freeze-out at 130 
MeV. We use an energy density weighted average for the 
computation of the eccentricity. For both V2 and e, we 
take into account the fact that the principle axes of the 
initial matter distribution are tilted with respect to the 
reaction plane. So we get non-zero values even for very 
central collisions, due to the random fluctuations. 

For all theoretical curves, the ratio V2/eis not constant, 
but increases substantially from peripheral towards cen- 
tral collisions - in agreement with the data. In our case, 
this increase is a core-corona effect: for peripheral colli- 
sions (small number of participating nucleons A'part); the 
relative importance of corona to core increases, and since 
the corona part does not provide any V2, one expects 
roughly |7J] 



!^ - f (M ^ ^ 

^ — ^ core V-' * part 7 ' ^ 



(31) 



with a monotonically increasing relative core weight 
/coro(-^Vpart), which varics between zero (very peripheral) 
and unity (very central). Comparing the theoretical 
curves in flg. [27l we see that most elliptical flow is pro- 
duced early, as seen by the dashed line, representing an 
early freeze out - at Tpo = 7h = 166 MeV. Adding fi- 
nal state hadronic rescattering leads to the full curve 
(full cascade) or the dotted one (only elastic scattering), 
adding some more 20 % to V2- The difference between 
the two rescattering scenarios is small, which means the 
effect is essentially due to elastic scatterings. Continuing 
the hydrodynamic expansion through the hadronic phase 
till freeze out at a low temperature (130MeV), instead of 
employing a hadronic cascade, we obtain a even higher 
elliptic flow, as shown by the thin line in fig. [571 and as 
discussed already in[l^ [l^ [75| . 

We now discuss the effect of the equation of state (see 
also [Z^). Using a (non-realistic) first-order equations of 
state (curve QIF from fig. [TT|) . one obtains considerably 
less elliptical fiow compared to the calculation using the 
the cross-over equation of state X3F, as seen in fig. 
Taking a wrong equation-of-state and a wrong treatment 
of the hadronic phase (thermally equilibrated rather than 
hadronic cascade) compensate each other, concerning the 
elliptical flow results. 

In our realistic (ideal) hydrodynamical treatment we 
get always an increase of the ratio of V2 over eccentricity, 
whereas it is also claimed that this variation is due to 
incomplete thermalization [77| . 

More detailed information is obtained by investigat- 
ing the (pseudo)rapidity dependence of the elliptical flow, 
for different centralities, as shown in flg. for Au-Au 
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Figure 28: Centrality dependence of the ratio of V2 over ec- 
centricity, for a full calculation, hydro & hadronic cascae, 
for a (non-realistic) first-order transition equations of state 
(dashed-dotted line) compared to the cross-over equations of 
state, the default case (full line, same as the one in fig. [27]). 
Points are data 1731. 
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Figure 29: Pseudorapidity distributions of the elliptical flow 
V2 for minimum bias events (upper left) and different central- 
ity classes, in Au-Au collisions at 200 GeV. Points are data 
[731, the different curves refer to the full calculation - hydro 
& cascade (full thick line), only elastic hadronic scatterings 
(dotted), no hadronic cascade at all (dashed), and hydrody- 
namic calculation till final freeze-out at 130 MeV (thin line). 



scattering at 200 GeV. Again we compare several scenar- 
ios: the full treatment, namely hydrodynamic evolution 
from flux tube initial conditions with early hadronization 
(at 166 MeV) and subsequent hadronic cascade, and the 
calculations with only elastic rescattering, or no hadron 
scattering at all. Also shown as thin line is the case 
where the hydrodynamic expansion is continued through 
the hadronic phase till freeze out at a low temperature 




(ISOMeV), instead of employing a hadronic cascade. The 
previously found observations are confirmed: at central 
rapidity, most flow develops early, the non-equilibrium 
hadronic phase gives only a moderate contribution. At 
large rapidities, however, the hadronic rescattering has a 
big relative effect on V2- Remarkable is the almost tri- 
angular shape of our V2 rapidity dependencies. This is 
partly due to the fact that the initial energy density is 
provided by flux tubes, each one covering a certain width 
in (space-time) rapidity, as indicated in fig. [26l A single 
elementary flux tube contributes a constant energy den- 
sity in a given interval, where the interval always contains 
rapidity zero. If (for a simple argument) the positive 
string endpoints were distributed uniformly in rapidity 
between zero and rjf'^^, the energy density would be of 
the triangular form 



(32) 



what we observe approximately. This initial shape in 
space-time rapidity rj^ seems to be mapped to the pseudo- 
rapidity dependence of V2 ■ 

Also important for this discussion is the fact that the 
relative corona contribution is larger at large rapidities 
compared to small ones. The corona contributes to par- 
ticle production (visible in rapidity spectra), but not to 
the elliptical flow. 

The above V2 results we obtained by averaging over 
transverse momenta pt , with the dominant contribution 
coming from small transverse momenta. The pt depen- 
dencies of V2 for different particle species is shown in fig. 
[501 (for minimum bias Au-Au collisions) and [3T] (for the 
20-60% most central Au-Au collisions), where we com- 
pare our simulations for pions, kaon, and protons with 
experimental data. We first look at the results for the 
transverse momentum dependence of V2 for the calcula- 
tions without hadronic cascade (w/o HC), i.e. the upper 
left plots in figs. [301 and [311 The pion and kaon curves 
are almost identical, the protons are shifted, due to an 
important corona contribution (considering only core, all 
three curves are on top of each other). Turning on the 
final state hadronic cascade (upper right plots) will pro- 
vide the mass splitting as observed in the data. Although 
this mass splitting was considered a great success of the 
hydro approach, it is in reality provided by the (non- 
thermal) hadronic rescattering procedure. It is this final 
state hadronic rescattering which is responsible for the 
fine structure of the pt dependence, although the magni- 
tude of the integrated V2 is produced in the early phase. 
The lower panel of the figs. [301 and [3T] shows a some- 
what different presentation of the same results: here we 
plot the scaled quantity V2/'nq versus the scaled kinetic 
energy [mt — rn)/nq^ where nq is the number of quarks of 
the corresponding hadron (2 for mesons, 3 for baryons). 
We show again the calculation without (left) and with 
(right) hadronic cascade. And surprisingly it is this final 
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Figure 30: The transverse momentum dependence of vi for 
pions (circles, full lines), kaons (squares, dashed lines), and 
protons (triangles, dotted lines) for minimum bias events in 
minimum bias Au-Au collisions at 200 GeV. The symbols refer 
to data(79l. [80|. the lines to our full calculations. 




Figure 31: The transverse momentum dependence of v-z for 
pions (circles, full lines), kaons (squares, dashed lines), and 
protons (triangles, dotted lines) for the 20-60% most central 
events in Au-Au collisions at 200 GeV. The symbols refer to 
datafsj, the lines to our full calculations. 



state hadronic rescattering which makes the three curves 
for pions, kaons, and protons coincide. At least in the 



small "Pt region considered here, the key for understand- 
ing "'i;2scaling" is the hadronic cascade, not the partonic 
phase. 



VII. GLAUBER OR COLOR GLASS INITIAL 
CONDITIONS 

There has been quite some discussion in the litera- 
ture concerning the possibility of increasing the elliptical 
flow when using Color Glass Condensate initial condi- 
tions rather then Glauber ones |82l. l83|. The latter ones 
are usually based on a simple Ansatz, assuming that the 
energy density is partly proportional to the participants 
and partly to the binary scatterings. 
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Figure 32: Initial energy density as a function of the radius r 
for azimuthal angles (j) = and (j> = '''"/2, from six randomly 
chosen flux tube initial conditions (full thin line: (j) = 0, dotted 
thin line: (p = n/2) and from Color Glass Condensate initial 
conditions (full line: = 0, dashed line: = 7r/2), for a 
semi-peripheral Au-Au collision. 

In our case, we compute partial cross sections, which 
gives us the number of strings (elementary flux tubes) 
per nucleon-nucleon collision. So we have as well contri- 
butions proportional to the binary scatterings (the string 
contributions), in addition to the remnant excitations, 
being proportional to the participants. On the other 
hand, we do consider high parton density effects, intro- 
ducing screening. In addition, the hydrodynamic expan- 
sion only concerns the core, and cutting off the corona 
pieces will produce sharper edges of the radial energy 
density distribution. In fig. I32[ we compare the energy 
density distributions as obtained from a CGC calcula- 
tion [3, [S^l , with six randomly chosen different events 
from our flux tube initial condition, after removing the 
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Figure 33: Initial energy density as a function of the radius r 
for azimuthal angles (j> = and (f> = n/2, from six randomly 
chosen flux tube initial conditions (full thin line: (j> = 0, dotted 
thin line: — n/2) and from Color Glass Condensate initial 
conditions (full line: (p — 0, dashed line: (j> = 7r/2), for a 
semi-peripheral Au-Au collision. 



corona. In fig. 1331 we compare the same distributions 
from the same same six individual events to calculations 
from Glauber initial conditions [isl . [s^l- Seeing these 
large event-by-event fluctuations, it is difficult to imag- 
ine that the differences between CGC results and Glauber 
are an issue when doing event- by- event treatment.. 



VIII. TRANSVERSE MOMENTUM SPECTRA 
AND YIELDS 

We have discussed so-far very interesting observables 
like two-particle correlations and elliptical flow. How- 
ever, we can only make reliable conclusions when we also 
reproduce elementary observables like simple transverse 
momentum (pt) spectra and the integrated particle yields, 
for identified hadrons. We will restrict the following pt 
spectra to values less than 1.5 GeV (2 GeV in some cases), 
mainly in order to limit the ordinate to three or at most 
four orders of magnitude, which allows still to see 10% 
differences between calculations and data. 

In the upper panel of fig. [Ml we show the pt spectra 
of 7r'''(left) and tt^ (right) in central Au-Au collisions, for 
rapidities (from top to bottom) of 0, 2, and 3. The middle 
panels show the transverse momentum / transverse mass 
spectra of tt"*" and w~ , for different centralities, and the 
lower panel the centrality dependence of the integrated 
particle yields per participant for charged particles and 
TT~ mesons. In fig. [3SJ we show the corresponding results 
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Figure 34: Production of pions in Au-Au collisions at 200 
GeV. Upper panel: transverse momentum spectra for central 
collisions at different rapidities (from top to bottom: 0, 2, 
3). The lower curves are scaled by factors of 1/2 and 1/4, 
for better visibility. Middle panels: transverse momentum 
(mass) distributions at rapidity zero for different centrality 
classes: from top to bottom: the 0-5%, the 20-30%, and the 
40-50% most central collisions. Lower panel: the centrality 
dependence of the integrated yields for charged particles and 
pions. The symbols refer to data [ssl - IS^ . the full lines to our 
full calculations, the dotted lines to the calculations without 
hadronic cascade. 



for kaons. In the upper panels, for the y = 2 and y = 3 
curves, we apply scaling factors of 1/2 and 1/4, for bet- 
ter visibility, all other curves are unsealed. We present 
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Figure 35: Same as fig. 1351 but for kaons. 
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always two calculations: the full one (full lines), namely 
hydrodynamic evolution plus final state hadronic cascade, 
and the calculation without cascade (dotted lines) . There 
is a slight increase of pion production in particular at low 
Pt during the hadronic rescattering phase, but the dif- 
ference between the two scenarios is not very big. We 
see almost no difference between between the calculation 
with and without hadronic rescattering in case of kaons. 
For both, pions and kaons, we observe a change of slope 
of the Pt distributions with rapidity. Concerning the cen- 
trality dependence, we observe an increase of the yields 
per participant. 

In fig. [551 a-nd 1371 we show pt spectra and central- 
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Figure 36: Production of lambdas (left) and antilambdas 
(right) in Au-Au collisions at 200 GeV. Upper panel: trans- 
verse momentum distributions at rapidity zero for different 
centrality classes: from top to bottom: the 0-5%, the 20- 
30%, and the 40-50% most central collisions. The lower curves 
are scaled by factors of 1/2, 1/4, and 1/8, for better visibil- 
ity. Lower panel: the centrality dependence of the integrated 
yields. The symbols refer to data (46| . the full lines to our 
full calculations, the dotted lines to the calculations without 
hadronic cascade. The thin line refers to a hydrodynamic cal- 
culation till final freeze-out at 130 MeV. 
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ity dependence of particle yields per participant, for the 
(multi) strange baryons A, A, S, and S. Same conven- 
tions as for the previous plots. Here we see a big effect 
due to rescattering: for the lambdas, the yields are not 
affected too much, but the pt spectra get much softer, 
when comparing the full calculation with the one with- 
out rescattering. Similarly the slopes for the S, and S 
get softer due to rescattering. 

We also show in the lower panels of figs. [36] and |37] 
the yields per participant in case of a hydrodynamic cal- 
culation till final freeze-out at 130 MeV (thin lines). We 
have almost no centrality dependence, in contrast to the 
significant increase seen in the data, for both, lambdas 
and xis. Such a full thermal scenario with late freeze-out 
is therefore incompatible with strange baryon data. 

For xis, the softening of pt spectra due to hadronic 
rescattering is more pronounced for the antiparticles ~ an 
absorption effect. Even the total integrated yields are af- 
fected: rescattering will reduce the S yields and increase 
the 5 yields with centrality. Maybe too much absorp- 
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Figure 38: Same as fig. 1371 but comparing the calculation 
without hadronic cascade (dotted) with the one witli only 
elastic hadronic rescattering (full thin line). 

tion? In fig. [38l we replace the full hadronic cascade by 
an option where only elastic rescattering is allowed (full 
lines). The dotted line refers to the calculation without 
rescattering, as in the previous plots. Here - by definition 
- the yields are unchanged, only the slopes are affected. 
It seems that this option reproduces the data better than 
the full cascade. 

In any case, the effect of rescattering decreases with 
decreasing centrality: the interaction volume simply gets 
smaller and smaller, reducing the possibility of rescatter- 
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Figure 39: Transverse momentum spectra of protons (left) 
and antiprotons (right) in Au-Au collisions at 200 GeV. Up- 
per panel: spectra for central collisions at different rapidities 
(from top to bottom: 0, 2, 3). The lower curves are scaled by 
factors of 1/2 and 1/4, for better visibility. Middle and lower 
panels: transverse momentum (mass) distributions at rapid- 
ity zero for different centrality classes: from top to bottom: 
the 0-5%, the 20-30%, and the 40-50% most central collisions. 
The symbols refer to datafi^. [ssl . [87| . the full lines to our 
full calculations, the dotted lines to the calculations without 
hadronic cascade. 

We finally discuss proton and antiproton production. 
When talking about spectra of identified hadrons, it is 
implicitly assumed that these spectra do not contain con- 
tamination from weak decays, so the experimental spec- 
tra should be feed-down corrected - which is not always 
the case. This is in particular important for protons, 
strongly affected by feed-down from lambda decays. So 
whenever we compare to data, we adopt the same defini- 
tions: in case of feed-down correction of the data, we sup- 
press weak resonance decays, and in case of no feed-down 
correction, we do let them decay. So for the following 
discussion, in case of the STAR data we compare to, pro- 
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tons are contrary to the pions not corrected, we include 
weak decay products. When comparing to PHENIX and 
BRAHMS data, we suppress weak decays. In fig. 1311 we 
show the the proton and antiproton transverse momen- 
tum spectra at different rapidities and different centrali- 
ties, for Au-Au collisions at 200 GeV. Again we show the 
full calculation (full lines) and the one without hadronic 
cascade (dotted lines). There is a huge difference be- 
tween the two calculations, so proton production is very 
strongly affected by the hadronic cascade. Not only the 
slopes change, also the total yields are affected. 

To summarize the above discussion on yields and pt 
spectra: an early hadronization at 166 MeV gives a rea- 
sonable description of the particle yields, which are not 
much affected by the hadronic final state rescattering, 
except for the protons. The main effect of the hadronic 
cascade is a softening of the pt spectra of the baryons. 



IX. FEMTOSCOPY 

All the observables discussed so-far are strongly af- 
fected by the space-time evolution of the system, nev- 
ertheless we investigate the momentum space, and con- 
clusions about space-time are indirect, as for example our 
conclusions about early hadronization based on particle 
yields and elliptical flow results. A direct insight into the 
space-time structure at hadronization is obtained from 
using femtoscopical methods (89l - l93| . where the study of 
two-particle correlations provides information about the 
source function S'(P,r'), being the probability of emit- 
ting a pair with total momentum P and relative distance 
r'. Under certain assumptions, the source function is re- 
lated to the measurable two-particle correlation function 
C(P,q) as 



^?(P,q) 



'5(P,r')|*(q',r')r 



(33) 



with q being the relative momentum, and where ^' is 
the outgoing two-particle wave function, with q' and r' 
being relative momentum and distance in the pair center- 
of-mass system. The source function S can be obtained 
from our simulations, concerning the pair wave function, 
we follow [oil, some details are given in appendix [FI 



As an application, we investigate tt^ 



correlations. 



Here, we only consider quantum statistics for ^ ^ no final 
state interactions, to compare with Coulomb corrected 
data. To compute the discretized correlation function 



C(Pi,qj), we do our event-by-event simulations. 



and compute for each event C-^ = Impairs I* (q'' ''O I ) 
where the sum extends over all 7r+ pairs with P and q 
within elementary momentum-space-volumes at respec- 
tively Pi and (\j . Then we compute the number of pairs 
Nij for the corresponding pairs from mixed events, be- 
ing used to obtain the properly normalized correlation 



function = C[j/Nij. The correlation function will be 
parametrized as 



C(P,q) = 

1 -I- A exp ( — i?out ^out ^ -^sidc 9sidc ^ -^i 



(34) 



2 

lono 



ftong) 5 



where "long" refers to the beam direction, "out" is par- 
allel to projection of P perpendicular to the beam, and 
"side" is the direction orthogonal to "long" and "out" 
[qsI - IqtI . In fig. HOI we show the results for the fit param- 
eters A, i?out) -Rsidc, and i?iong, for five different central- 
ity classes and for four fc^ intervals defined as (in MeV): 
KT1= [150,250], KT2= [250,350], KT3= [350,450], 
KT4= [450, 600], where kr of the pair is defined as 



= ^ (IPT(pionl) 



PT(pion2)|) . 



(35) 



Despite what appears in j98| , this is the correct defini- 
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Figure 40: Femtoscopic radii i?out, i?sido, and i?iong, as well 
as A as a function of rriT for different centralities (0-5% most 
central, 5-10% most central, and so on). The full lines are the 
full calculations (including hadronic cascade), the stars data 

M 



tion of kr used by STAR in their analysis [9; 
suits are plotted as a function of rriT = \/k^ 



The re- 
^. The 



model describes well the radii, the experimental lambda 
values are sightly below the calculations, maybe due to 
particle misidentification. Both data and theory provide 
lambda values well below unity, maybe due to pions from 
long-lived resonances. Concerning the to^ dependence 
of the radii, we observe the same trend as seen in the 
data [9II: all radii decrease with increasing tot, and the 
radii decrease as well with decreasing centrality. This 
can be traced back to the source functions, shown in fig 
l¥T] These source functions are by definition the distribu- 
tions of the distances a;i(pionl) — a;j;(pion2) of the pairs. 
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Figure 41: The source functions as obtained from our sim- 
ulations, for three different centraUties (0-5% most central, 
10-20% most central, and 30-50% most central), representing 
the distribution of the space separation of the emission points 
of the pairs, in the "out" - "side" - "long" coordinate system, 
in the longitudinal comoving frame. The different curves per 
plot correspond to the different kr bins, see text. 




x(fm) 

Figure 42: The mean transverse momentum component px of 
tt"'" as a function of the x coordinate of the emission point. 
Also shown is the number of produced 7r+ as a function of x. 
The different curves refer to different centralities: 0-5% = full 
line, 10-20% = dashed, 30-50% = dotted. 




where Xi are coordinates of the emission points. We use 
the "out" - "side" - "long" coordinate system, and the 
longitudinal comoving reference frame. To account for 
the fact that only small values of the magnitude of the 
relative momentum |q| provide a non-trivial correlation, 
we only count pairs with |q| < 75 MeV. The different 
curves per plot correspond to the different values of kx 
bins: the upper curve (full red) correspond to KTl, the 
second curve from the top (dashed blue) correspond to 
KT2, and so on. In other words, the curves get narrower 
with increasing kx, which is perfectly consistent with the 
decreasing radii in fig. |40l Concerning the centrality de- 
pendence, the curves get narrower with decreasing cen- 
trality, in agreement with decrease of radii with decreas- 
ing centrality seen in fig. UHl 

The reason for the decrease of radii with mx is the 
strong space-momentum correlation. In fig. I42[ we show 
the average of produced tt"*" mesons as a function of 
the X coordinate of their formation positions, for different 
centralities. Clearly visible is the strong x ~ px correla- 
tion, being typical for radial flow. Also visible in the 
figure is the smaller spatial extension for peripheral com- 
pared to central collisions. To illustrate this phenomenon, 
we show in fig. 221 a situation of completely radial trans- 
verse momentum vectors, who's magnitudes increase with 
increasing distance from the center. We consider two 
pairs of momentum vectors, a and b at some distance 
ri as well as c and d at some distance r2 <ri. We have 
chosen the pairs such that the magnitude of their differ- 



Figure 43: Radial flow effect on mt dependence of femtoscopic 
radii. 



ences is the same (and "small"), to mimic the fact that 
only pairs with small relative momentum are relevant for 
the HBT analysis. The spatial distance between the two 
momentum vectors c and d is bigger than the one for the 
pair a and b. due to the fact that the latter vectors are 
longer than the former ones (|a| ~ \h\ > \c\ sa |d|). In this 
way we understand the connection between increasing nit 
and decreasing space separation. 

We now consider two other scenarios: the calculation 
without hadronic cascade (final freeze out at 166 MeV), 
and the fully thermal scenario, where we continue the hy- 
drodynamical evolution till a late freeze-out at 130 MeV 
(and no cascade afterwards either). In figs. l44l and |45| 
we see a similar space-momentum correlation as for the 
complete calculation in fig. 1421 the mean transverse mo- 
mentum components p^ is roughly a linear function of the 
transverse coordinate x, in the region where the particle 
density is non-zero. The maximum mean px is smaller 
in the no-cascade case, and bigger in the fully thermal 
case, as compared to the complete calculation. Interest- 
ing are the dn/dx distributions: the no-cascade results 
(with early hadronization) are much narrower than the 
full thermal ones. The complete calculation of fig. |42] 
is in-between, in the sense that the plateau of the dn/dx 
distribution is similar to the no-cascade case, but the tails 
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Figure 44: Same as fig. 1421 but for tlie calculation without 
hadronic cascade. 
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Figure 45: Same as fig. 1421 but for the full thermal scenario 
(freeze-out at 130 MeV. 



widths of the single particle source functions dn/dx. For 
large values of rout) the "full thermal" scenario and the 
one "without cascade" coincide - as do the shapes of the 
tails of the single particle source functions. A similar be- 
havior is found for all the source functions, as shown in 
figs. ST] and 051 where we plot the source functions for 
the "full thermal" and the "without cascade" scenarios. 
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Figure 47: Same as fig. 1411 but for a calculation without 
hadronic cascade. 
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Figure 46: The source function for the "out" coordinate for 
the three scenarios: complete calculation, with hadronic cas- 
cade (full line), calculation without hadronic cascade and 
therefore final hadronization at 166 MeV (dashed), and full 
thermal scenario with hydrodynamic evolution till the final 
freeze-out at 130 MeV (dotted). 

In fig. we compare the source functions for the 
three scenarios, namely the complete calculation, the cal- 
culation without hadronic cascade, and the full thermal 
scenario with hydrodynamic evolution till the final freeze- 
out. For small values of Tout, the "complete calculation" 
and the "full thermal" one coincide - as do the total 
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Figure 48: Same as fig. 1411 but for a calculation where the 
hydro evolution is continued till freeze-out at 130 MeV (being 
the final freeze-out, no cascade afterwards). 

The above discussion is important to understand the 
results concerning the femtoscopic radii for the different 
scenarios. The fitting procedure used to obtain the fem- 
toscopic radii is based on the hypothesis that the source 
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functions are Gaussians, the fit is therefore blind con- 
cerning the non-Gaussian tails. Due to the fact that the 
source function from the complete calculations and the 
full thermal scenario are identical apart from the tails, 
we expect similar results for these two scenarios, whereas 
the calculation without cascade should give smaller radii. 
This is exactly what we observe in fig. 031 where we show 
femtoscopic radii for the calculations without hadronic 
cascade (full line) and with hydrodynamical evolution till 
final freeze-out at 130 MeV (dashed). We observe always 
a decrease of the radii with mT^ but the dependence is 
somewhat weaker as compared to the data. But the mag- 
nitude in case of "no cascade" is very low compared to 
the two other scenarios, which are relatively close to each 
other, and to the data. Here the radii do not allow to 
discriminate between two scenarios which have neverthe- 
less quite different source functions. This is a well-known 
problem, and ther e are met hods to go beyond Gaussian 
parameterizations but we will not discuss this 

any further. 
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Figure 49: Same as fig. 1401 but the calculations are done 
without hadronic cascade (full line) or with a hydrodynamic 
evolution through the hadronic phase with freeze-out at 130 
MeV (dashed). 

Although the Gaussian parameterizations represent 
only an incomplete information about the source func- 
tions, the centrality and transverse momentum depen- 
dence of the radii is nevertheless very useful. It is a 
necessary requirement for all models of soft physics to 
describe these radii correctly. There has been for many 
years an inconsistency, referred to as "HBT puzzle" (6^ . 
Although hydrodynamics descibes very successfully ellip- 
tical flow and to some extent particle spectra, one can- 
not get the femtoscopic radii correctly, when one uses 
"simple" hydrodynamics. Using transport models (and 
an event-by-event treatment) may help [H]. In [6a|. it 
has been shown that the puzzle can actually be solved 



by adding pre-equilibrium flow, taking a realistic equa- 
tion of state, adding viscosity, using a more compact or 
more Gaussian initial energy density profile, and treat- 
ing the two-pion wave fun ction more accuratly. It has 
also been shown [l06l - [l08t that using a Gaussian initial 
energy density profile, an early starting time (equivalent 
to initial flow), and a cross-over equation of state, and a 
late sudden freeze-out (at 145 MeV) helps to descibe the 
femtoscopic radii, and to so me extent the spectra. 

The scenario in [l06l - IT08t is compatible with our sce- 
nario "hydrodynamical evolution till final freeze-out at 
ISOMeV", which allows us to get the femtoscopic radii 
correctly (see fig. H^ . as well as some V2 results and 
some spectra. One cannot describe, however, yields and 
spectra of lambdas and xis. 



X. SUMMARY AND CONCLUSIONS 

We presented a realistic treatment of the hydro- 
dynamic evolution of ultrarelativistic heavy ion colli- 
sions, based on flux-tube initial conditions, event-by- 
event treatment, use of an efficient (3+l)D hydro code 
including flavor conservation, employment of a realistic 
equation-of-state, use of a complete hadron resonance ta- 
ble, and a hadronic cascade procedure after an hadroniza- 
tion from thermal matter at an early time. 

Such an approach is able to describe simultaneously 
different soft observables such as femtoscopic radii, par- 
ticle yields, spectra, and V2 results. One obtains in a 
natural way a ridge structure when investigating ArjAcj) 
correlations, without adding a particular mechanism. 

Considering such a multitude of observables, a clear 
picture of the collision dynamics emerges: a hydrody- 
namic evolution starting from initial flux-tube structures, 
till hadronization at an early time in the cross-over region 
of the phase transition, with subsequent hadronic rescat- 
terings being quite important to understand the shapes 
of particle spectra. 
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Appendix A: Pomeron structure 

We define a so-called profile function function G asso- 
ciated to a Pomeron exchange as 



G{b) = ^2Imf (6), 



(Al) 



with T being the Fourier transform of the Pomeron ex- 
change scattering amplitude T, 

m = 4^ I d?q^ e-*«-^>(0, (A2) 

using t = —q\. 

There are two contributions, a soft and a semi-hard 
one. The energy-momentum dependence of the semi-hard 
profile function may be expressed in terms of light cone 
momentum fractions as 

G'sGmi(^PEi ^Pe) — -^part (^pe) -^part (^Pe) ^(^PE^Pe); 

(A3) 

where the vertex function impart is given as 



using 



impart (a;) = a-px'^" , 



"F = S^^^^7/M /3f = ec - Opa 



(A4) 



(A5) 



with parameters ec, 7?i, apart, and with 

w(4e^pe) - / dx+dx^ J dtY,E\Ml,x+)E'{Ml,'. 



da 



The variables x^ are light cone momentum fractions. The 
QCD evolution function is computed in the usual way 
based on the DGLAP equations, 



dFP"' 



d\uQ^ 



E 

k 



z 2tt 



with the initial condition 



E'^Zb {Q' = Ql,x)=SJ^6{l-x). 



(A8) 



(A9) 



Here P^{z) are the usual Altarelli-Parisi splitting func- 
tions. One introduces the concept of "resolvable" parton 
emission, i.e. an emission of a final (s-channel) parton 
with a finite share of the parent parton light cone mo- 
mentum (1 — z) > e = p5_res /Q^ (with finite relative 
transverse momentum p']_ = Q^{1 ~ z) > p\j.^g ) and use 
the so-called Sudakov form factor, corresponding to the 
contribution of any number of virtual and unresolvable 
emissions (i.e. emissions with (1 — z) < e) , 




dq^ 
,2 



dz^Pl:iz) 



A^Q2,Q2)=exp. 

iQl r 

(AlO) 

This can also be interpreted as the probability of no re- 
solvable emission between Q§ and . Then -Eqcd can 

be expressed via -Eq™p) , corresponding to the sum of any 
number (but at least one) resolvable emissions, allowed 
by the kinematics: 

+ &^c^{QlQ\x), (All) 
where i?Q™p) (QQ,(3^,a;) satisfies the integral equation 
EZcT,{QlQ^x) (A12) 



Q2 



dQl 



E 



dz as 
~2tt 



/ + — + — i\ 
(.ajpgXpgXgXg s,t). 



dt 



(A6) 



The indices i and j refer to parton fiavors. Alp is the fac- 
torization scale (here Mp = tu/s). The quantity datj/dt 
is the hard Born parton-parton scattering cross section, 
and E'^{A1p,xe) the so-called complete evolution func- 
tion, being a convolution of the soft and the QCD evolu- 
tion, 



+ A^{QlQl)^Pp{x)\ A'"(Q?,Q2). 

Here Pj^{z) are the Altarelli-Parisi splitting functions for 
real emissions, i.e. without (5-function and regularization 
terms at z 1. Eq. (|A12p can be solved iteratively, see 
I25i]. 

We define the soft contribution Gsoft (•?,&) as |25fl 



G soft {s,b) 



with 



27f 



part 



Asoft(s/so) 



exp 



4Asoft(s/so) 
(A13) 



i?XM|,4) = ^ f dxt,,dx' 

7. -J 



± 

QCD 



(A7) 



Kotiz) = 2i?pa,.t + a soft Inz, 



(A14) 



Eaoiti^toft) -^QCd(^^I'j2;qcd)'^(2^E ~ ^^ft^QCD) 



with parameters agoft, a'soft 
So = IGeV^. 



7part, -Rpart: ^nd a scale 
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Appendix B: Solving hydrodynamic equations 

The algorithm is based on the Godunov method: one 
introduces finite cells and computes fluxes between cells 
using the (approximate) Riemann problem solution for 
each cell boundary. A relativistic HLLE solver is used to 
solve the Riemann problem. To achieve more accuracy 
in time, a predictor-corrector scheme is used for the sec- 
ond order of accuracy in time, i.e. the numerical error is 
Oidt^), instead of 0{dt'^). To achieve more accuracy m 
space, namely a second order scheme, the linear distribu- 
tions of quantities (conservative variables) inside cells are 
used. The conservative quantities are {e+p*v'^)/ {1 — v'^), 
(e +p) * v/{l — v"^) . 

We rewrite equations in hyperbolic coordinates. These 
coordinates are suitable for the dynamical description at 
ultrarelativistic energies. It is convenient to write the 
equations in conservative form, the conservative variables 
axe 



/{Qr 



Q = 



1 Qr\ 




7^(6 +p)Wx 


Qx 




Qy 






Qn 






Qb 




inB 


Qs 




ins 


\QqJ 


\ inq J 



(Bl) 



where ub, ns, nq are the densities of the conserved quan- 
tities _B, S, and Q. The components Qm are conservative 
variables in the sense that the integral (discrete sum over 
all cells) of Qm gives the total energy, momentum, and 
the total B, S, and Q, which are conserved up to the 
fluxes at the grid boundaries. The velocities in these ex- 
pressions are defined in the "Bjorken frame" related to 
velocities in laboratory frame as 



..lab 



,,lab 



coshy 



cosh(y - rjs) 
coshy 



^ ^ cosh(j/ - rjs) 
Vr, = tanh(2; - r]s) 



(B2) 



where y = + vl^^)/{l - vl"^^)] is the longitudinal 

rapidity of the fluid element, r/s = ^ \n[{t + z) / {t — z)] is 
space-time rapidity. The full hydrodynamical equations 
are then 



dr 



/ Qr\ 




( Qr 






/ V(p • v) \ 


Qx 




Qx 






dxP 


Qy 




Qy 






dyP 


Qrj 


+v • 


Qr] 




v + 




Qb 




Qb 









Qs 




Qs 









\Qq) 




\Qq 


J 




V J 


quantities 






fluxes 



(B3) 



-p)(l + 

Qx/t 

Qy/r 

2Qr,/T 

Qb/t 
Qs/r 

Qq/t 



sources 



with V = {dx, dy, ^d^). 

We base our calculations on the finite- volume approach 
: we discretize the system on a fixed grid in the calcula- 
tional frame and interpret Q^, j^j. as average value over 
some space interval AVijk, which is called a cell. The 
index n refers to the discretized time. 

The values of QJ^ j^j, are then updated after each time- 
step according to the fiuxes on the cell interface during 
the time-step At„. One has the following update formula 



At 



At 



AX2 

At 
A^ 



Axi 

(-^«(,i+i/2),fc 



iP{i+l/2),jk + F(i-l/2),jk) 



F, 



) (B4) 



{Ftj.(k+l/2) + -Fij,(fc-l/2)), 



where F is the average fiux over the cell boundary, the 
indexes -1-1/2 and —1/2 correspond to the right and the 
left cell boundary in each direction. This is the base of 
the Godunov method |l09| . which also implies that the 
distributions of variables inside a cell are piecewise linear 
(or piecewise parabolic etc, depending on the order of 
the numerical scheme), which forms a Riemann problem 
at each cell interface. Then the flux through each cell 
interface depends only on the solution of a single Riemann 
problem, supposing that the waves from the neighboring 
discontinuities do not intersect. The latter is satisfied 
with the Gourant-Friedrichs-Lewy (GFL) condition |llOl |. 

To solve the Riemann problems at ea ch cell interface, 
we use the relativistic HLLE solver which approxi- 

mates the wave profile in the Riemann problem by a single 
intermediate state between two shock waves propagating 
away from the initial discontinuity. Together with the 
shock wave velocity estimate, in this approximation one 
can obtain an analytical dependence of the flux on the 
initial conditions for the Riemann problem, which makes 
the algorithm explicit. 

We proceed then to construct a higher-order numerical 
scheme: 

• in time: the predictor- corrector scheme is used for 
the second order accuracy in time, i.e. the numeri- 
cal error is 0{dt^), instead of 0{dt^) 

• in space: in the same way, to achieve the second 
order scheme, the linear distributions of quantities 
(conservative variables) inside cells are used. 
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Some final remarks: 

At each time-step, we compute and sum the fluxes 
for each cell with all its neighbors and update the value 
of conservative variables with the total flux. Thus, we 
do not use operator splitting (dimensional splitting) and 
thus avoid the numerical artifacts introduced by this 
method, e.g. artificial spatial asymmetry. 

To treat grid boundaries, we use the method of ghost 
cells. We include 2 additional cells on either end of grid 
in each direction, and set the quantities in these cells 
at the beginning of each time-step. For simplicity, we 
set the quantities in ghost cells to be equal to these in 
the nearest "real" cell, thus implementing non-reflecting 
boundary conditions (outflow boundary). This physically 
correspond to boundary which does not reflect any wave, 
which is consistent with expansion into vacuum. 

In our simulations we deal with spatially finite systems 
expanding into vacuum. Thus the computational grid 
in Eulerian algorithm must initially contain both system 
and surrounding vacuum. To account for the finite ve- 
locity of the expansion into the vacuum, which equals c 
for an infinitesimal slice of matter on the boundary, we 
introduce additional (fioating-point) variables in each cell 
which keep the extent of matter expansion within a cell, 
having the value unity for the complete cell, zero for a 
cell with vacuum only. The matter is allowed to expand 
in the next vacuum cell only if the current cell is filled 
with matter. 



with 

Hi = B,HB + QilJ-Q + Si^j,s, (C6) 

where hb, fJ-s^ fJ-Q a-re the chemical potentials associ- 
ated to B, 5, Q, and Bi, Si, Qi are the baryon charge, 
strangeness, and the electric charge of i-th hadron state, 
gi = (2 Ji -I- 1) is degeneracy factor. 

For large baryon chemical potential the EoS correction 
for the deviations from ideal gas due to particle interac- 
tions becomes more important. We employ this correc- 
tion in a form of an excluded volume effect, like a Van 
der Waals hard core correction. According to this pre- 
scription, 

i 

fl^^H^-Vi-p. (C8) 

If one supposes equal volume Vi = v for all particle 
species, then the correction can be computed as a solution 
p(r, fiB, fJ-Q, fJ-s) of a fairly simple, however transcenden- 
tal equation, 

p{T, tiB,t^Q,l^s) = /°'*^(T, fiB,f^Q, f,s)e-^P^^'^-'^'^-'^^y^ 

(C9) 

We take the value v w 1.44 fm^, which corresponds to 
the hard core radius r = 0.7 fm. 



Appendix C: resonance gas 



Appendix D: Ideal QGP 



Whereas for hadronization we employ the correct quan- 
tum statistics, we use the Boltzmann approximation for 
the calculation of the equation of state. This is reason- 
able even for pions at zero chemical potential, the ex- 
cluded volume correction at nonzero chemical potentials 
is considerably bigger than the difference coming from 
quantum statistical treatment. We account for all well 
known hadrons made from u, d, s quarks from the PDG 
table For energy density, pressure and net charges we get 



^m^y + _Ak,{^) cxp(M./r) 



P = 



T 



T 



riB = 



^^m2T2.X2(f )-exp(M,/T) 

i 

^i?,^m?T.if2(^)-exp(^,/T) 



ns 



= E0'&"^?^-^^2(^)-exp(M,:/T) 



'27r2 



T 
T 



(CI) 
(C2) 

(C3) 

(C4) 

(C5) 



In this ideal phase, matter is made from massless m, d 
quarks and massive s-quark (+antiquarks). Due to the 
possibility of a large strange quark chemical potential, 
comparable to its mass mg = 120 MeV which is taken 
in our calculations, we perform the integration of the 
strange quark contribution to thermodynamic quantities 
exactly, without Boltzmann or zero-mass approximation. 
So we have 



91 

67r2 
9i 

67r2 



4 



60 



(Dl) 



60 

2 



+ PsiT, p.s) +ps{T,Hs) + 
with pj(T, lis) = ps{T, -Hs), and 



90 



^ rj-i f-OQ 

27r^ Jo 



1- 



-exp 



Ms 

T 



dp, 



(D2) 

where we use the degeneracy factors gi — 6 for light 
quarks, gg = 16 for gluons, and a bag constant B = 
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0.38 GeV/fm^. Quark chemical potentials are 



A'i'^i - Pi we get 



3 67r2 



1 91 



[2^.3 + 27^'^luT' ^VrfT'] - (D8) 

(D9) 



ns = - [7is(r,/xs) - ?ij(r, -//s)] 

with e5(r, fis) = es(r, -^s), and 



oo 


oo 



cxp ( ^V?^^ 



1 



expl^^j^T 



s T 

- I + 1 



7,2 - ii^ 



-dp, 
(DIO) 
dp. 
(Dll) 



Appendix E: Plasma hadronization 

We parametrize the hadronization hyper-surface = 
x^{t, ip, if) as 

a;" = TC0sh?7, = rcosyj, x^ = rsin(/5, x^ = rsinhyy, 

(El) 

with r = r{T,ip,r]) being some function of the three pa- 
rameters T, if, rj. The hypersurface element is 



dx" dx"" (9x^ , , , 
= ^A.-'^^ Qt dri (^2) 



with e'^'^'^^ = —£fiuKX = 1- Computing the partial deriva- 
tives dx^ /da, with a = r, (p, 77, one gets 





fD3) 


dSo = 




(D4) 


dSi = 


MS- 


(D5) 


dS2 = 




= Ts + 


= 


:+B 


(D6) 


Cooper-Fry( 




(D7) 





Or Or 1 

-r—r cosh 77 + r— — sinhyy > drdipdri, (E3) 
or 1777 J 

dr 1 

— — r sin <j9 + r r cos (j9 > drdLpdri, (E4) 
J 

— — — T cos 09 -f- r r sin 09 > drdipdri, (E5) 
o</3 J 

dr , dr 1 

?'— rsinh?? — r — cosh?/ > drd(j9(i?7.(E6) 

OT OTj J 



d^p 



d^fiP^fiup), 



with M being the flow four-velocity in the global frame, 
which can be expressed in terms of the four- velocity u in 
the "Bjorken frame" as 



u = u cosh rj + u sinh 77 , 



u " sinh 77 + M cosh 7/ 



r, 3, 



(E7) 
(E8) 
(E9) 
(ElO) 



In a similar way one may express p in terms of p in the 
Bjorken frame. Using 7 = {t " and the flow velocity = 
u^'/j, we get 



dn 



dyd(j)dp±^ 



~ r 

rrp 



(Ell) 



with p '' = p^ cos i^H-p2 sin and p* = p^ sin ip — p'^ cos 
being the radial and the tangential transverse momentum 
components. Our Monte Carlo generation procedure is 
based on based on the invariant volume element moving 
through the FO surface, 



with 



dV* — dY^ii^u'^ = w drdifidr], 



dr dr ^ dr , 

-r—T + rrv + -r— ru — r—v" 
dr dtp drj 



(E12) 



(E13) 



and with = cosLp + simp and — sin 93 — 
cos If being the radial and the tangential transverse 
flow. Freeze out is the done as follows (equivalent to 
Cooper- Frye): the proposal of isotropic particles produc- 
tion in the local rest frame as 



diM^ad^p* dV*ME*), 



(E14) 
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is accepted with probability 



adV*E* 



(E15) 



In case of acceptance, the momenta are boosted to the 
global frame. 



Appendix F: Pair wave function for femtoscopy 
applications 



In case of identical particles, we use 



1 



$(q',r') = -^(0(k',r')±</>(-k',r')), (Fl) 



and for non-identical particles 

$(q',r')=<^(k',r'), 



(F2) 



with k' = q'/2. In the simplest case, neglecting final 
state ineteractions, one has simply 



-k',r') =cxp(-ik',r'). 



(F3) 



otherwise the non-symmetrized wavefunction is given as 
(see eq. (89) of H) 



-k',r') = eMiSc) VAJji) 



(F4) 



exp(-*k', r') F{-^^^, 1, tO + Uk')^^ 



with £, = k'r' + q'r' , p = k'r', r] = (k'a)-^. The 
quantity a ~ (//ziZ2e^)^^is the Bohr radius of the 
pair, in case of pion-pion one has 387 fm. Furthermore, 
5c = argr(l + irj) is the Coulomb s-wave phase shift, 
Ac{ri) = 27r?7 (exp(27r77) — 1) ^ is the Coulomb penetra- 
tion factor, 

F(a, 1, z) = 1 + az/V."^ + a{a + \)z'^ 12"? + ... (F5) 

is the confluent hypergeometric function, 

G(p,77) =P(p,77)+277pB(p,77) (F6) 
X [ln|277p| + 2C-l + x(?7)], 



with the Euler constant C =^ 0.5772, and 

oo oo 

B{p,Ti)=Y,Bs. P{p.il) = Y,Ps. (F7) 

with Bq = 1, Bi= rjp, Pq = 1, Pi = 0, and 

(n + l)(n + 2)B„+i = 2r]pB„ - p^B^-u (F8) 
n{n + l)P„+i = 27?pF„ - p^P,,^i (F9) 
- (2?! + l)277pB„. 

The function x is given as 

x{v)^h{Tj)+tAc{v)/{2v), (FIO) 

where h is expressed in terms of the digamma function 
ip{z) = r'{z)/r{z) as 

Hv) = \ [V^("?) + VX-^^y) - Hv")] ■ (Fii) 

The amplitude fc can be written as 

Uk') = /(fc')Mc(r?), (F12) 

where /(fc') is the amplitude of the low energy s-wave 
elastic scattering due to the short range interaction renor- 
malized by the long-range Coulomb forces. We may write 



/c(fc') = K 



,^-1 2x(ry) 



(F13) 



with [ni 



_ 2l_ ^th - So ^ ^ / 2fc' ^ 



s - So 



/Sth 



(F14) 



s - \/mf + fc'2 j , sth = (mi + m2)2, (F15) 

with the parameters as given in [ll2j . 
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